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Abstract 

We present lattice results for the ground state energy of a spin-1/2 fermion system in the unitary 
limit, where the effective range of the interaction is zero and the scattering length is infinite. We 
compute the ground state energy for a system of 6, 10, 14, 18, and 22 particles, with equal numbers 
of up and down spins in a periodic cube. We estimate that in the limit of large number of particles, 
the ground state energy is 0.25(3) times the ground state energy of the free Fermi system. 
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I. INTRODUCTION 



The unitary limit of spin- 1/2 or two- component fermions has attracted much recent 
experimental and theoretical interest from several branches of physics. The term unitary 
limit or unitary regime refers to the idealized limit where the effective range of the interaction 
is zero and the scattering length is infinite. Much of the interest has been spurred by 
experimental advances in the trapping of ultracold atomic Fermi gases. Starting with a 
dilute Fermi gas, where the effective range of the interaction is negligible compared with 
the interparticle spacing, one can reach the unitary limit by tuning the scattering length to 
infinity using a Feshbach resonance 

In nuclear physics the unitary limit is directly relevant to the properties of cold dilute 
neutron matter. The neutron scattering length is roughly —18 fm while the effective range is 
2.8 fm. Therefore the unitary limit is approximately realized when the interparticle spacing 
is about 5 — 10 fm, roughly 0.5% — 5% of normal nuclear matter density, which is believed 
to be relevant to the physics of the inner crust of neutron stars jjj. 

At zero temperature in the unitary limit there are no dimensionful parameters other than 
the particle density. Therefore the ground state energy of the system should obey a simple 
relation E° = £E 0,free , where E 0, is the ground state energy of a non-interacting Fermi gas 
and £ is a dimensionless constant. Recent experiments have measured the expansion of 6 Li 
in the unitary limit released from a harmonic trap (f| 7L Based on a Thomas- Fermi model, 
the measured values for £ are 0.51(4) |J and 0.321}° 3|. The discrepancy between these 



two recent measurements and with larger values for £ reported in earlier experiments 
suggest further experimental work is needed, as well as a better theoretical understanding 
of Thomas- Fermi theory and other local density approximations in the unitary limit. 

There have been several recent analytic calculations of £ using various techniques such as 
BCS saddle point and variational approximations, Pade approximations, mean field theory 
with pairing, and dimensional expansions 0, Q, Q, Q, Q|. The values for £ vary from 
roughly 0.3 to 0.6. Fixed-node Green's function Monte Carlo calculations have found £ to 

Qn 
and 0.42(1) [15j. A recent estimate based on Kohn-Sham theory for the two- 

fermion system in a harmonic trap yields a value of 0.42 Q|. Recently there have also been 

simulations of spin-1/2 fermions on the lattice in the unitary limit . 
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19, 201. While 



the simulations are at nonzero temperature, the results can be extrapolated to provide an 
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estimate for £. The results of |3| produced a value for £ similar to the fixed-node results, 
while Q established a bound, 0.07 < £ < 0.42. 

In this paper we describe a new calculation of £ on the lattice at zero temperature and 
fixed particle number. In contrast with Green's function Monte Carlo, the method we 
describe is free from the fermion sign problem. Because of this we can eliminate systematic 
errors due to fermion nodal constraints. Using a new algorithm which combines endpoint 
correlation function importance sampling and non-local Monte Carlo updating, we measure 
the ground state energy for 6, 10, 14, 18, and 22 particles, with equal numbers of up and 
down spins in a periodic cube. We estimate that in the limit of large number of particles, 
the ground state energy is 0.25(3) times the ground state energy of the free Fermi system. 



II. LATTICE FORMALISM 

We refer to the state with N spin-up fermions and N sp in-down fermions as an N, N 
state. We use the same lattice action as described in Since we will be working 

at fixed particle number we can set the chemical potential to zero. Throughout we use 
dimensionless parameters and operators, which correspond with physical values multiplied 
by the appropriate power of the spatial lattice spacing a. For quantities in physical units 
we use the superscript 'phys'. The lattice action has the form 

n,i 

- h J2 [c*(n) Ci (n + l s ) + c*(n) Ci (n - /,)] + ~ ^ ( s ( H )f ■ W 

n,l B ,i ft 

Here n labels the sites of a 3 + 1 dimensional lattice, l s (s = 1, 2, 3) are the spatial lattice 
unit vectors, is a temporal lattice unit vector, and i labels the two spin components of the 
fermion. The temporal lattice spacing is a t , and at = a t /a is the ratio of the temporal to 
spatial lattice spacing. We have also defined h = a t /(2m), where m is the fermion mass in 
lattice units. 

Our choices for the physical values of the fermion mass and lattice spacings are irrelevant 
to the universal physics of the unitary limit. Nevertheless it is convenient to assign some 
concrete values to these parameters. The values we choose are motived by the dilute 
neutron system. We use a fermion mass of 939 MeV and lattice spacings a = (50 MeV) -1 , 
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at = (24 MeV) -1 . Our spatial geometry is a periodic cube of length L lattice units. The 
Grassmann fields are denoted by Ci(n). s(n) is an auxiliary Hubbard-Stratonovich field 
which upon integration reproduces the attractive contact interaction between up and down 
spins. The interaction coefficient for these lattice spacings is computed by summing two- 



particle scattering bubble diagrams. The details of the calculation are given in [2jJ . Taking 
the scattering length a scat t — > oo, we find in the unitary limit, C phys = —1.203 x 10~ 4 MeV -2 . 
In order to compute the ground state energy, we consider the correlation function 



J N,N 



(t) = (*%, N \e- m \K,N), (2) 



where the initial/final state l^jviv) is the normalized Slater determinant ground state for 
the free N, N particle system. We refer to t as the Euclidean time. We define 

E N>N (t) = -^[\nZ N , N (t)\. (3) 

Then as t — ► +oo, E^^it) converges to E^ N , the ground state energy of the interacting 
N, N particle system. The only assumption is that the ground state has a nonvanishing 
overlap with the ground state of the non-interacting system. 

For the non-interacting system we find E^ e ^(t) = E^]^, where E^^ is the free Fermi 
ground state energy. It is useful to define the dimensionless function 

E NjN (t) 

^N,N \t) - o.free • W 

In the unitary limit ^Ar,jv(t) depends only on the dimensionless combination —j2- We can 
define the Fermi energy, 



E F = ^ = — (3n 2 ) 2/3 ~ 7 596^ (5) 
2m 2m mL 2 ' 



where kp is the Fermi momentum. In the unitary limit we can regard £jv,iv(0 as a function 
of Ept, and for sufficiently large N it should approach a common asymptotic form. 

The conversion of the Grassmann lattice action to a worldline formalism at fixed particle 
number has been detailed in (2^. We use the same transfer matrix derived there, except in 
this case we keep the auxiliary Hubbard-Stratonovich field and calculate the integral over 
auxiliary field configurations, 

Z N<N {t) cxjDs e-^(*)f (*% iN \ Te- H ^ \^ N ) . (6) 



The advantage of the auxiliary field formalism is that H(s) consists of only single-body 
operators interacting with the background auxiliary field. We can therefore compute the 
full N, iV-body matrix element as the square of the determinant of the single-particle matrix 
elements, 

(*%,N\ Te ~ H{s)t \*%) « [detM(,,t)] 2 , (7) 
M ij (s,t) = (p i \Te- H M\ Pj ), (8) 

where i,j go from 1 to N. The states \pj) are single-particle momentum states comprising 
our Slater determinant initial/final state. The square of the determinant arises from the 
fact that the single-body operators in H(s) are the same for both up and down spins. Since 
the square of the determinant is nonnegative, there is no sign problem in this formalism. 

Our notation, Te~ H ^ 1 , is shorthand for the time-ordered product of single-body transfer 
matrices at each time step, 

Te -H(s)t = M ( Lt _!).___. M ( nt ) . _ _ _ . M ( X ) . M ( ), (9) 

where L t is the total number of lattice time steps and t = L t a t . If the particle stays at the 
same spatial lattice site from time step n t to n t + 1, then the corresponding matrix element 
of M{n t ) is 

e V=CET t s(il)+22L^ 1 _ 6h y ^ 

If the particle hops to a neighboring lattice site from time step n t to n t + 1 then the corre- 
sponding matrix element of M(n t ) is h. All other elements of M(n t ) are zero. 



III. ENDPOINT IMPORTANCE SAMPLING AND HYBRID MONTE CARLO 

While there is no sign problem, the calculation of Z^^{t) is not straightforward. Due to 
the large coupling strength in the unitary limit, fluctuations in [det M(s, t)] 2 are very large. 
We deal with this problem by sampling configurations according to the weight 

6X P j"^ E ( S ^)) 2 + 21 °§ H det M (Mend)|] j , (11) 

where t en d is the largest Euclidean time at which we wish to measure Z N;N (t). It is most 
efficient to sample configurations using a non-local updating method called hybrid Monte 
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Carlo |23[. A description of the hybrid Monte Carlo method as applied to the grand 
canonical ensemble of spin-1/2 neutrons can be found in j^]. For this fixed particle number 
simulation we compute molecular dynamics trajectories for 



(12) 



where p(n) is the conjugate momentum for s(n) and 

v(s) = lJ2 ( s ^)) 2 ~ 2 lo s [' det M ( s ' tend ) 



(13) 



The steps of the algorithm follows. 



Step 1: Select an arbitrary initial configuration s°( 



n 



Step 2: Select p°(n) according to the Gaussian random distribution 



Step 3: Let 



P(jr(n)) oc exp 



p°(n) = p°(n) — 



dV(s) 



ds(n) 



(14) 



(15) 



for some small positive e. In computing the derivative of V, we use the fact that 



dV(s) 
ds(n) 



2 ^ ddet M dM H 
8 ^ n ' ~ detM ^ dM kl ds(n) 



s(n)-2j2[M- 1 } 

k,i 



dM kl 
lk ds(n) 



(16) 



Step 4: For j = 0, 1, N - 1, let 



s J+ (n) = s 3 (n) + spin), 



Step 5: Let 



p l+ (n) = p 3 (n) — e 



p N {n)=p N (n) + 



dV(s) 
ds{n) 



dV(s) 
ds{n) 



(17) 
(18) 

(19) 
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Step 6: Select a random number r G [0, 1). If 

r < exp [-H(s N , p N ) + H(s°, p )] (20) 

then let 

s°(n) = s N (n). (21) 
Otherwise leave s° as is. In either case go back to Step 2. 

For each configuration the observables that we calculate are 

[det M(s,t))> (22) 
[detM( S ,t end )] 2 ' 

for t < t cn d. By taking the ensemble average of 0(s,t) we are able to calculate 

Z N,N(t) ^ 
Zn,N (tend) 

Fluctuations in the observable 0(s, t) are manageable in size so long as t en d — t, the temporal 
separation from the endpoint, is not too large. This is typically not a problem since only 
a small window in t near t en d is needed to calculate — [In Z^^{t)] at t = t cn< j. For each 
simulation we compute roughly 2 x 10 5 hybrid Monte Carlo trajectories, split across four 
processors running completely independent trajectories. Averages and errors are computed 
by comparing the results of each processor. 



IV. ERROR ESTIMATES 



We use double precision arithmetic to compute det M(s, t e nd) and 0(s,t). We carefully 
monitor any systematic errors produced by double precision round off error and exceptional 
configurations. To do this we introduce a small positive parameter e and reject any hybrid 
Monte Carlo trajectories which generate a configuration with 

IdetM^e* JJ |M«| . (24) 

i=l,...,N 

We then take the limit e — > + to determine if poorly condition matrices make any detectable 
contribution to our observable. We consider values for e as small as 1CT 6 . If as we take 
e — > + any systematic error can be detected above the stochastic error level, then we throw 
out the measurement and do not include it in the final results. The error bars we present 
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are therefore estimates of the total error for each lattice system. There are no additional 
errors other than the lattice spacing dependence. 

Errors due to finite lattice spacing can be estimated by the spread of £jv,iv(^) at fixed 
N for different lattice lengths L. By fitting the data for different L to a single function 
£,N,N{t), we can compute the uncertainty in the fit parameters and estimate the error of the 
measurement in the continuum limit. Therefore the error bars on our final results include 
both systematic and stochastic errors for each lattice system as well as systematic errors 
due to lattice spacing dependence. 



V. RESULTS 



As a first test we consider the 1, 1 particle system. In this case the free Fermi ground 
state energy E^ ee is zero. So instead of dividing by E^'f 66 , we measure the dimensionless 
quantity mL 2 Ei^{t). In Table 1 mL 2 Ei^(t) is shown for various t and L. We also show 
the results for the dimensionless interacting ground state energy, mL 2 E\ x , as computed by 
summing lattice bubble diagrams. 

Table 1: mL 2 Ei^(t) for various t and L 



L 


I2a t 


24a* 


36a^ 


48a t 


60a t 


mL 2 E\ x 


4 


-2.7(1) 


-3.4(1) 


-3.4(1) 


-3.4(1) 


-3.5(1) 


-3.57 


5 


-2.5(1) 


-3.1(1) 


-3.5(1) 


-3.5(1) 


-3.5(1) 


-3.60 


6 


-2.0(1) 


-2.7(1) 


-3.1(1) 


-3.6(1) 


-3.5(1) 


-3.62 



For large t we see that mL 2 Ei^(t) matches mL 2 E\ l within error bars. For large L the 
lattice discretization error should be negligible and mL 2 E\ x should approach the unitary 
limit result, 

4-K 2 d x ~ 4tt 2 (-0.095901) ~ -3.7860, (25) 
where d\ is the large scattering length expansion coefficient in a periodic cube as defined in 



24j . We see that mL 2 E\ x is within a few percent of the unitary limit value for L = 4,5, 6. 



We have performed hybrid Monte Carlo simulations of the N = 3, 5, 7, 9 systems with 
lattice lengths L = 4, 5,6 and the N = 11 system with lattice lengths L = 5,6. In Figs. 
EJ 121 El HI Q we show £,N,N(t) as a function of Ept for iV = 3, 5, 7, 9, 11 respectively. 
In all cases we find agreement in £jv,jv(£) for different values of L, indicating the physics 
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FIG. 1: f 3)3 (t) versus E F t for L = 4,5,6. 
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FIG. 2: f 5)5 (t) versus E F t for L = 4,5,6. 



signature of the unitary limit. The agreement for different values of L shows that to a good 
approximation L is the only dimensionful scale in the system. The lack of any significant 
corrections due to varying L indicates that the lattice discretization error is small and that 
the scattering length is essentially infinite. This provides a non-trivial check within the 
context of the many body system that the physics is consistent with the unitary limit. We 
find that detuning the coupling constant away from the unitary limit value by as little as 
5% in either direction produces a measurable disagreement in £jv,jv(i) for different values of 
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FIG. 3: f 7)7 (t) versus Ep* for L = 4,5,6. 
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FIG. 4: f 9)9 (t) versus £ F t for L = 4,5,6. 



L. We discuss this further in the next section. 

We fit £,N,N{t) at large Ept to the functional form b + cexp(— 5 • Ept) to determine the 
asymptotic value as E F t — > oo. This is the asymptotic form one expects for very large E F t, 
where 5 ■ E F is the gap in the energy spectrum above the ground state energy E^ N . If 
there is a non-negligible coupling to gapless phonon modes, then 5 should go to zero as we 
increase the number of particles and volume. However it is not clear which excited states 
have a non-negligible overlap with our free particle initial/final state. It is also unclear if 
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FIG. 5: £ii,ii(i) versus Ept for L = 5,6. 



we are probing at sufficiently large Ept and with sufficient accuracy to determine 5 reliably. 
We hope to resolve these questions with future studies. For our purposes here we consider 
S only as a fit parameter to determine the asymptotic value for £,N,N{t) at large Ept. 

The ratio of the ground state energy of the interacting system to the free particle ground 
state energy is given by the large Ept limit of £jv,jv(i)- The results for these ratios are shown 
in Table 2. 

Table 2: Results for E% N f E^ r ^ e 



N = 3 


N = 5 


N=7 


N = 9 


N = 11 


0.19(2) 


0.24(2) 


0.28(2) 


0.23(2) 


0.25(2) 



From Table 2 we see that the energy ratio for the smallest system, N = 3, is somewhat 
lower than the rest. However the ratios for N > 5 are close to a central value of about 0.25. 
Assuming no large changes to this ratio for iV > 11, we estimate that £ = 0.25(3). 

In Fig. IHlwe show simultaneously all of the plots of £jv,at(0 as a function of E F t. We 
see again that the smallest system, N = 3, is somewhat different from the rest, falling off 
with a faster exponential at large Ept. The curves for N = 5,7, 9, 11 seem quite similar in 
both shape and asymptotic value, suggesting that the curves are already close to the large 
N limit. 
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FIG. 6: £,N,N(t) versus -Epi for iV = 3,5,7,9, 11. We have drawn lines showing the central value 
and error bounds for the estimate £ = 0.25(3). 



VI. DISCUSSION AND ADDITIONAL CROSS-CHECKS 

Our result £ — 0.25(3) disagrees with the fixed-node Green's function Monte Carlo results 
£ = 0.44(1) [l4j and 0.42(1) 15j. However it is consistent with the requirement that the 
fixed-node value sets a variational upper bound on the energy. One explanation of the 
disagreement could be the difficulty in determining a lower bound for £ using Green's function 
Monte Carlo. This requires removing the nodal constraint and overcoming a sign problem 
that scales exponentially with N. But in any case the disagreement of our final result and 
published fixed-node results suggests further tests of the robustness of our method. We 
consider three additional cross-checks in this section. 

The first test we consider is whether or not the non-quadratic lattice dispersion relation 
is significantly affecting our results. While deviations were found in nonzero temperature 
simulations [19], the fact that %N,N(t) is the same for different values of L suggests that this 
is not the case at zero temperature. However if we are measuring true continuum limit 
behavior, it should be possible to reproduce the same results using a different lattice action 
with the same continuum limit. To carry out this cross-check we consider an 0(a 2 )-improved 
action with next-to-nearest neighbor hopping that reproduces the continuum dispersion 
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FIG. 7: 0(a 2 )-improved and standard lattice action results for £3,3(2) versus Ept for L = 4,5,6. 

relation p 2 /(2m) up to 0(p 6 ). Instead of the standard action (JTJ, we use the improved 
action 



c*(n) Ci (n + 6) 



2 



- \h ^ Jc* (n)cj(n + f s ) + c*(n)Ci(n - l s ) 



n,LA 



12 



h J2 \<( H ) C ^ H + + c*(n)ci(n - 2l s ) ( s (™)) 



(26) 



As before we calculate the interaction coefficient by summing two-particle scattering bubble 
diagrams. We find in the unitary limit, C phys = —1.581 x 10~ 4 MeV -2 . 

The results for N = 3 are shown in Fig. We have taken the standard lattice action 
results shown in Fig. [Hand superimposed the new 0(a 2 )-improved lattice data. There is 
almost no difference between the standard and 0(a 2 )-improved lattice results. This suggests 
that we are measuring continuum limit behavior. 

The second test we consider is whether the interaction coefficient can be retuned to recover 
the fixed-node value for £ while approximately retaining the physics of the unitary limit. 
In order to test this possibility we consider the N = 5 system. We tune the interaction 
coefficient so that for L = 4 we get lim^oo £5,5(2) = 0.42. By trial and error we find that 
this occurs for C phys = —1.08 x 10 -4 MeV -2 , which is about 90% of the unitary limit value. 
In the two-particle system this corresponds with a scattering length of a scatt — —7.0 fm 
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FIG. 8: £5,5 (t) versus E F t for L = 4, 5, 6 and C phys = -1.08 x 1(T 4 MeV" 2 . 
~ -0.035 MeV -1 . 

We now simulate the system for L = 5 and L = 6. The results are shown in Fig. |SJ 
The plots for £5,5(2) do not agree for different L, and we find limt_ +DO £5,5(2) = 0.42, 0.50, 
0.57 for L = 4, 5,6 respectively. From these results it appears that limt-+oo £5,5(2) = 0.42 
is inconsistent with the unitary limit. In fact we can use this data to estimate the finite 
scattering length correction. Using a scatt ~ —0.035 MeV -1 we find 

lim £5,5(2) « 0.24 - 0.6 • k~ F l a; c l tt . (27) 

A more thorough study of the finite scattering length correction to £ will be presented in 
future work. 

The third test we consider is whether the same ground state energy can be extracted 
using different initial/final states. We test this using the N = 7 system with L — 4. 
We define ^77) as the state we get by taking ^77) and removing the four fermions at 
momentum (p x ,P y ,Pz) = (^O^fr) an d setting them at higher momentum, (p x ,P y ,Pz) = 
(iff, iff, 0). We define |tf£ >7 ) in the same way, except we put only two fermions at the 
higher momentum while keeping the total momentum equal to zero. 
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FIG. 9: Comparison of £7,7 (i), £7,7 (t), and £7,7 (*)■ 



We define the following quantities using the excited Slater determinant states: 

ZUt) = (* 7>7 | e"« |^ 7 ) , Zl 7 (t) = <* 2 7J | e~ m |* 2 7j7 ) , (28) 
4r(*) = ~J t MJ, r (f)] > *?,r(*) = "| M,r(*)] - ( 2 9) 

£7,7 (*) = pO.frce ' ^7,7^) = „Q,free ' ( 30 ) 
^7,7 ^7,7 

The comparison of £7,7(2), £7 7 (t), and £77(2) is shown in Fig. The convergence of all three 
results at large Ept suggests that we are measuring the true ground state energy at large 
E F t. 



VII. SUMMARY 

We have measured the ground state energy for N, N spin- 1/2 fermions in the unitary 
limit in a periodic cube. Our results at N = 3, 5, 7, 9, 11 suggest that for large N the 
ratio of the ground state energy to that of a free Fermi gas is 0.25(3). Our result lies in 
the middle of the bound 0.07 < £ < 0.42 given in j^J. Although it is inconsistent with 
the fixed-node Green's function Monte Carlo results £ = 0.44(1) Q and 0.42(1) Q, it is 
consistent with the requirement that the fixed-node value sets a variational upper bound 
on the energy. We have cross-checked the robustness of our results by comparing with an 
0( a 2 ) -improved action, detuning away from the unitary limit to measure deviations, and 

15 



using different initial/final states to reproduce the same ground state energy. In this paper 
we have not addressed the questions of the existence of superfluidity, long range order, or 
the pairing gap. However the methods presented here shows some promise at probing some 
of these interesting questions. 
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